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S ! Abstract 

The reheating of the universe after hybrid inflation proceeds through the nucleation 
■ and subsequent collision of large concentrations of energy density in the form of 

' bubble-like structures moving at relativistic speeds. This generates a significant frac- 

tion of energy in the form of a stochastic background of gravitational waves, whose 
time evolution is determined by the successive stages of reheating: First, tachyonic 
preheating makes the amplitude of gravity waves grow exponentially fast. Second, 
I . bubble collisions add a new burst of gravitational radiation. Third, turbulent mo- 

^jq' tions finally sets the end of gravitational waves production. From then on, these 

waves propagate unimpeded to us. We find that the fraction of energy density today 
in these primordial gravitational waves could be significant for GUT-scale models of 
inflation, although well beyond the frequency range sensitivity of gravitational wave 
observatories like LIGO, LISA or BBO. However, low-scale models could still produce 
a detectable signal at frequencies accessible to BBO or DECIGO. For comparison, 
we have also computed the analogous background from some chaotic inflation models 
and obtained similar results to those of other groups. The discovery of such a back- 
ground would open a new observational window into the very early universe, where 
the details of the process of reheating, i.e. the Big Bang, could be explored. Thus, it 
could also serve as a new experimental tool for testing the Inflationary Paradigm. 
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1 Introduction 

" I 

■ Gravitational waves (GW) are ripples in space-time that travel at the speed of light, and whose emission 

^ I by relativistic bodies represents a robust prediction of General Relativity. Theoretically, it is expected 

that the present universe should be permeated by a diffuse background of GW of either an astrophysical or 
cosmological origin [T]. Fortunately, these backgrounds have very different spectral signatures that might, 
in the future, allow gravitational wave observatories like LIGO [5^, LISA BBO [|] or DECIGO [F, to 
disentangle their origin Jj. Unfortunately, the weakness of gravity will make this task extremely difficult, 
requiring a very high accuracy in order to distinguish one background from another. 

There are, indeed, a series of constraints on some of these backgrounds, coming from the anisotropics 
in the Cosmic Microwave Background (CMB) 6 , from Big Bang nucleosynthesis 17! or from millisecond 
pulsar timing [Sj. Most of these constraints come at very low frequencies, from 10~^^ Hz to 10^^ Hz, while 
present and future GW detectors (will) work at frequencies of order 10~^ — lO'^ Hz. If early universe first 
order phase transitions [HI HH] or cosmic turbulence jTT] occurred around the electro-weak (EW) scale, 
GW detectors could have a chance to measure the corresponding associated backgrounds. However, if 
those processes occurred at the GUT scale, their corresponding backgrounds will go undetected by the 
actual detectors, since these cannot reach the required sensitivity in the high frequency range of lO'' — 10^ 
Hz. There are however recent proposals to cover this range |12[ I13j . which may become competitive in 
the not so far future. 
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Cosmological observations seem to suggest that something hke Inflation must have occurred in the very 
early universe. Approximately scale-invariant density perturbations, sourced by quantum fluctuations 
during inflation, seem to be the most satisfying explanation for the CMB anisotropics. Together with 
such scalar perturbations one also expects tensor perturbations (GW) to be produced, with an almost 
scale- free power spectrum |14| . The detection of such a background is crucial for early universe cosmology 
because it would help to determine the absolute energy scale of inflation, a quantity that for the moment 
is still uncertain, and would open the exploration of physics at very high energies. 

However, in the early universe, after inflation, other GWB could have been produced at shorter 
wavelengths, in a more 'classical' manner rather than sourced by quantum fluctuations. In particular, 
whenever there are large and fast moving inhomogeneities in a matter distribution, one expects the emis- 
sion of GW. At large distances from a source, the amplitude of the GW is given by hij ~ GQij/c'^r, with 
Qij the quadrupole moment of the mass distribution. The larger the velocity of the matter distribution, 
the larger the amplitude of the radiation produced. However, because of the weakness of gravity, in order 
to produce a significant amount of gravitational radiation, it is required a very relativistic motion (and 
high density contrasts) in the matter distribution of a source. Fortunatelly, this is indeed believed to be 
the situation at the end of inflaton, during the conversion of the huge energy density driving inflation 
into radiation and matter, at the so-called reheating of the Universe |15j . i.e. at the Big Bang. 

Note that any background of GW coming from the early universe, if generated below Planck scale, 
immediately decoupled upon production and, whatever their spectral signatures, they will retain their 
shape throughout the expansion of the Universe. Thus, the characteristic frequency and shape of the 
GWB generated at a given time should contain information about the very early state of the Universe 
in which it was produced. Actually, it is conceivable that, in the not so far future, the detection of these 
GW backgrounds could be the only way we may have to infer the physical conditions of the Universe at 
such high energy scales. However, the same reason that makes GW ideal probes of the early universe — 
the weakness of gravity — is responsible for the extreme difficulties we have for their detection on Earth. 

In Refs. [inilTT] we described the stochastic background predicted to arise from reheating after hybrid 
inflation. Here we will review the various processes involved in the production of such a background. In 
the future, this background could serve as a new tool to discriminate among different inflationary mod- 
els, since reheating in each model would give rise to a different GWB with very characteristic spectral 
features. The details of the dynamics of preheating depend very much on the model and are often very 
complicated because of the non-linear, non-perturbative and out-of-equilibrium character of the process 
itself. However, all the cases have in common that only specific resonance bands of the fields suffer an 
exponential instability, which makes their occupation numbers grow by many orders of magnitude. The 
shape and size of the spectral bands depend very much on the inflationary model. If one translates 
this picture into position-space, the highly populated modes correspond to large time-dependent inho- 
mogeneities in the matter distributions which acts, in fact, as a powerful source of GW. For example, in 
single fleld chaotic inflation models, the coherent oscillations of the inflaton during preheating generates, 
via parametric resonance, a population of highly occupied modes that behave like waves of matter. They 
collide among themselves and their scattering leads to homogenization and local thermal equilibrium. 
These collisions occur in a highly relativistic and very asymmetric way, being responsible for the genera- 
tion of a stochastic GWB [H Ull [Ml HD 112] with a typical frequency today of the order of lO"^ - 10^ Hz, 
corresponding to the present size of the causal horizon at the end of high-scale inflation. 

However, there are models like hybrid inflation in which the end of inflation is sudden [23j and the 
conversion into radiation occurs almost instantaneously. Indeed, hybrid models preheat very violently, 
via the spinodal instability of the symmetry breaking field that triggers the end of inflation, irrespective 
of the couplings that this field may have to the rest of matter. Such a process is known as tachyonic 
preheating [24. i25j and could be responsible for copious production of dark matter particles [26j , lepto 
and baryogenesis [27], topological defects [24], primordial magnetic fields [28], etc. In Ref. [25], it was 
shown that the process of symmetry breaking in hybrid preheating, proceeds via the nucleation of dense 
bubble-like structures moving at relativistic speeds, which collide and break up into smaller structures 
(see Figs. 7 and 8 of Ref. [H]). We conjectured at that time that such collisions would be a very strong 
source of GW, analogous to the GW production associated with strongly first order phase transitions [1] . 
As we will show here, this is indeed the case during the nucleation, collision and subsequent rescattering of 
the initial bubble-like structures produced after hybrid inflation. During the different stages of reheating 
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in this model, gravity waves are generated and amplified until the Universe finally thermalizes and enters 
into the radiation era of the Standard Model of Cosmology. From that moment until now, this cosmic 
GWB will be redshifted as a radiation-like fluid, totally decoupled from any other energy-matter content 
of the universe, such that today's ratio of energy stored in these GW to that in radiation, could range 
from U^^h"^ ~ 10~®, peacked around / ~ 10^ Hz for the high-scale models, to O^^/i^ ~ 10~^^, peacked 
around / ~ 1 Hz for the low-scale models. 

Finally, since the first paper by Khlebnikov and Tkachev [TS], studing the GWB produced at reheat- 
ing after chaotic infiation, there has been some developments. The idea was soon extended to hybrid 
inflation in Ref. [Hj. It was also revisited very recently in Ref. [501 [H] for the X(j)'^ and m^0^ chaotic 
scenarios, and reanalysed again for hybrid inflation in Refs. [16[|17| . using the new formalism of tachyonic 
preheating [211 US]. Because of the increase in computer power of the last few years, we are now able to 
perform precise simulations of the reheating process in a reasonable time scale. Moreover, understanding 
of reheating has improved, while gravitational waves detectors are beginning to attain the aimed sensitiv- 
ity [2j. Furthermore, since these cosmic GWBs could serve as a deep probe into the very early universe, 
we should characterize in the most detailed way the information that we will be able to extract from 
them. 



2 Gravitational Wave Production 

Our main purpose here is to study the details of the stochastic GWB produced during the reheating 
stage after hybrid inflation (sections 2 and 3). Nevertheless, we also study more briefly the analogous 
background from reheating in some chaotic models (section 4). Thus, in this section we derive a general 
formalism for extracting the GW power spectrum in any scenario of reheating within the (flat) Friedman- 
Robertson- Walker (FRW) universe. The formalism will be simplified when applied to scenarios in which 
we can neglect the expansion of the universe, like in the case of Hybrid models. 

A theory with an inflaton scalar field x interacting with other Bose fields 4'a, can be described by 

>C = Id^xd^'x + \d,.^ad^ctya + - 1^(0, X) (1) 

with R the Ricci scalar. For hybrid models, we consider a generic symmetry breaking 'Higgs' field $, 
with Nc real components. We can take = ^ 0^ = |0p/2, with a running for the number of Higgs' 
components, e.g. — 1 for a real scalar Higgs, Nc — 2 for a complex scalar Higgs or A^^, 4 for a SU{2) 
Higgs, etc. The effective potential then becomes 

F(0, x) = ^ - v'f + g\M' + \l^'x' ■ (2) 
For chaotic scenarios, we consider a massless scalar field (f> interacting with the infiaton x via 

v{xA)^\g\^4>'' + v{x), (3) 

with V{x) the infiaton's potential. Concerning the simulations we show in this paper, we concentrate in 
the TVc = 4 case for the hybrid model and consider a potential V{x) = ^X^ for chaotic scenario. 
The classical equations of motion of the inflaton and the other Bose fields are 

X + 3Hx^-V\+^^0, 0^ + 3i/</)a - -V20, + — = (4) 

a-^ ox a o<Pa 

with H ~ a I a. On the other hand, GW are represented here by a transverse-traceless (TT) gauge- 
invariant metric perturbation, /ly, on top of the flat FRW space ds^ = — rft^ + a^it) {Sij + hij) dx^dx^ , 
with o(t) the scale factor and the tensor perturbations verifying dihij = hu = 0. Then, the Einstein 
field equations can be splitted into the background and the perturbed equations. The former describe 
the evolution of the flat FRW universe through 

^ x' + ^{^x)' + ^l + ^{^M' (5) 



AttG 3a'^' 3a2 

•iH^ .o 1 , .9 1 



r + ^(vxr + + ^(v</>a)' + 2i-(x,0) (6) 
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where any term in the r.h.s. of (O and (O, should be understood as spatially averaged. 

On the other hand, the perturbed Einstein equations describe the evolution of the tensor perturba- 
tions |35l as 



% + 3Hh,j - = 16nGIl^j , (7) 

with 9^11^ = Ha = 0. The source of the GW, 11^ , contributed by both the inflaton and the other scalar 
fields, will be just the transverse-traceless part of the (spatial-spatial) components of the total anisotropic 
stress-tensor 

T^^ = [d^xdi^X + d^,(j)adu4>a + gt,i>{.C - (p))] /a^, (8) 

where >C(x, (pa) is the lagrangian ([T]) and {p) is the background homogeneous pressure. As we will explain 
in the next subsection, when extracting the TT part of ([H]), the term proportional to g^i, in the r.h.s 
of ([5]), will be dropped out from the GW equations of motion. Thus, the effective source of the GW will 
be just given by the TT part of the gradient terms S^tX^i/X + d^'f'adiy't'a- 

2.1 The Trans verse- Traceless Gauge 

A generic (spatial-spatial) metric perturbation Shij has six independent degrees of freedom, whose contri- 
butions can be split into [35j scalar, vector and tensor metric perturbations Shij = ip Sij+E^ij -\-F(^ij^+hij , 
with diFi = and dihij = ha = 0. By choosing a transverse-traceless stress-tensor source Hij, we can 
eliminate all the degrees of freedom (d.o.f.) but the pure TT part, hij, which represent the only physical 
d.o.f which propagate and carry energy out of the source (i.e. GW). Thus, taking the TT part of the 
anisotropic stress-tensor, we ensure that we only source the physical d.o.f. that represent GW. 
Let us switch to Fourier space. The GW equations ([7]) then read 

h^j{t,k)+3Hh,j{t,k) + —h,jit,k) = 167rGn,y(t,k), (9) 

where k = |k|. Assuming no GW at the beginnig of reheating (i.e. the end of inflation t^), the initial 
conditions are hij{te) = hij{te) = 0, so the solution to Eq. ([9]) for t > will be just given by a causal 
convolution with an appropriate Green's function G(t,t'), 



/iy (t,k) = ISttG^ dt' G{t,t')n^j{t',k) . 



(10) 



Therefore, all we need to know for computing the GW is the TT part of the stress-tensor, Hij, and the 
Green's function G{t' , t). However, as we will demonstrate shortly, we have used a numerical method by 
which we don't even need to know the actual form of G{t',t). To see this, let us extract the TT part of 
the total stress-tensor. Given the symmetric anisotropic stress-tensor T^i^ ([S]), we can easily obtain the 
TT part of its spatial components in momentum space, ny (k). Using the spatial projection operators 
Pi] = ^ij - f^i^ji with ki = ki/k, then ^ (k) Ay,;„i(k)T;„,.(k), where 

A,,-i™(k) = (Pa(y)PjM - ip„(k)Pz„(k)) . (11) 

Thus, one can easily see that, at any time t, kiHij (k, t) — n^(k, t) — 0, as required, thanks to the identities 

Pij kj — and Pij Pjm — Pim . 

Note that the solution (fTO|) is just linear of the non-traceless nor-transverse tensor Ty ([8]). Therefore, 
we can write the TT tensor perturbations (i.e. the GW) as 

h,j{t, k) = A,JJ„^{k)ul„^{t, k), (12) 

with Uij {t, k) the Fourier transform of the solution of the following equation 

ilij + SHiiij — ^V^Uij — WnGTij . (13) 
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This Eq. (fT^ is nothing but Eq. ([7]), sourced with the complete Tij ([5]), instead of with its TT part, 
Uij. Of course, Eq. contains unphysical (gauge) d.o.f.; however, in order to obtain the real physical 
TT d.o.f. hij, we can evolve Eg. p^ in configuration space, Fourier transform its solution and apply the 
projector (|lip as in p2p . This way we can obtain in momentum space, at any moment of the evolution, 
the physical TT d.o.f. that represent GW, hij. Whenever needed, we can Fourier transform back to 
configuration space and obtain the spatial distribution of the gravitational waves. 

Moreover, since the second term of the r.h.s of the total stress-tensor T^j is proportional to gij = 
6ij + hij, see ([S]), when aplying the TT projector (fTTj) . the part with the Sij just drops out, simply 
because it is a pure trace, while the other part contributes with a term — (£ — {p))hij in the l.h.s of 
Eq.®. However, (£ — (p)) is of the same order as the metric perturbation ^ 0{h), so this extra term 
is second order in the gravitational coupling and it can be neglected in the GW Eqs. This way, the 
effective source in Eq. is just the gradient terms of both the inflaton and the other scalar fields, 

= (V.xVjX + V,<?!)aVj(/)a)/a2. (14) 

Therefore, the effective source of the physical GW, will be just the TT part of ([H]), as we had already 
mentioned before. 

We have found the commuting procedure proposed (i.e. the fact that we first solve Eq. (jl3p and 
secondly we apply the TT projector to the solution (jl2p . and not the other way around), very useful. 
We are able to extract the spectra or the spatial distribution of the GW at any desired time, saving a 
great amount of computing time since we don't have to be Fourier transforming the source at each time 
step. Most importantly, with this procedure we can take into account backrcaction simultaneously with 
the fields evolution. 

In summary, for solving the dynamics of reheating of a particular inflationary model, we evolve 
Eqs. dH) in the lattice, together with Eqs. ©-([I]), while for the GWs we solve Eq. (fT3)) . Then, only 
when required, we Fourier transform the solution of Eq. (|13p and then apply ((T2)) in order to recover the 
physical transverse-traceless d.o.f representing the GW. From there, one can easily build the GW spectra 
or take a snapshot of spatial distribution of the gravitational waves. 



2.2 The energy density in GW 

The energy-momentum tensor of the GW is given by |36j 

where hij are the TT tensor perturbations solution of Eq. ([7]). The expectation value (...)y is taken 
over a region of sufficiently large volume V = to encompass enough physical curvature to have a 
gauge-invariant measure of the GW energy-momentum tensor. 
The GW energy density will be just = ton, so 
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J d='xAy(t,x)/i„(t,x) = J d^kk^{t,k)h*j{t,k), (16) 



where in the last step we Fourier transformed each hij and used the definition of the Dirac delta. We can 
always write the scalar product in p6p in terms of the (Fourier transformed) solution u/,„ of the Eq. (jl3p , 
by just using the fact that Ajj_;„iA;„i = ^ij.rs- This way, we can express the GW energy density as 

^•"^ " 32ttGL^ / ^^^'^ / ^^0M^i^)^^jit,k)Uhn{t,'^)■ (17) 



From here, we can also compute the power spectrum per logarithmic frequency interval in GW, normalized 
to the critical density pc, as fi^w ~ I ^ ^gw(/) i where 

- ^ I '^"A..™(k)..(i,k).L(i,k) (18) 

We have checked explicitely in the simulations that the argument of the angular integral of IjlSp is 
independent of the directions in k-space. Thus, whenever we plot the GW spectrum of any model, we 
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will be showing the amplitude of the spectrum (per each mode k) as obtained after avaraging over all the 
directions in momentum space, 

" 8G^ (Ay,i™(k)^y(i,k)^L(i,k))^^ (19) 
with EE 3L/ /da 

Finally, we must address the fact that the frequency range, for a GWB produced in the early universe, 
will be redshifted today. We should calculate the characteristic physical wavenumber of the present 
GW spectrum, which is redshifted from any time t during GW production. So let us distinguish four 
characteristic times: the end of inflation, te] the time when GW production stops; the time tr when 
the universe finally reheats and enters into the radiation era; and today, ^o- Thus, today's frequency /o 
is related to the physical wavenumber kt at any time t of GW production, via /o = (ai/ao)(fcf/27r), with 
ao and at, the scale factor today and at the time t, respectively. Thermal equilibrium was established 
at some temperature Tr, at time tr > t. The Hubble rate at that time was MpHr — (87r/3)/9r, with 
Pr — grTT^Tr /30 the relativistic energy density and gr the effective number of relativistic degrees of 
freedom at temperature Tr- Since then, the scale factor has increased as Gr/ao — {gQ,s/gr.sY^'^{To/Tr), 
with gi^s the effective entropic degrees of freedom at time ti, and To today's CMB temperature. Putting 
all together, 

where we have used the fact that the physical wave number kf at any time t during GW production, is 
related to the comoving wavenumber k through kt = {ae/a,t)k with the normalization = 1. 

From now on, we will be concerned with hybrid inflation, leaving chaotic inflation for section 4. Within 
the hybrid scenario, we will analyse the dependence of the shape and amplitude of the produced GWB 
on the scale of hybrid inflation, and more specifically on the v.e.v. of the Higgs field triggering the end 
of inflation. Given the natural frequency at hand in hybrid models, m — y/Xv, whose inverse sets 
the characteristic time scale during the first stages of reheating, it happens that as long as w ^ Mp, the 
Hubble rate H ^ \/X{v^ / Mp) is much smaller than such a frequency, H <ti m. Indeed, all the initial 
vacuum energy pQ gets typically converted into radiation in less than a Hubble time, in just a few m^^ 
time steps. Therefore, we should be able to ignore the dilution due to the expansion of the universe during 
the production of GW, at least during the first stages of reheating. Our approach will be to ignore the 
expansion of the Universe, such that we fix the scale factor to one, a — 1. As we will see later, neglecting 
the expansion of the Universe for the time of GW production, will be completely justified a posteriori. 

The system of equations that we have to solve numerically in a lattice for the hybrid model are 

iiy - V^Uy = 167rGTy (21) 

X-V\+(.g'|0p + ^2)x = O (22) 

4>a - V^4>a + (.g'x' + A|0P - m^) 0, = (23) 

with Tij given by Eq. p^ with the scale factor a — 1. We have explicitly checked in our computer 
simulations that the backreaction of the gravity waves into the dynamics of both the inflaton and the 
Higgs fields is negligible and can be safely ignored. We thus omit the backreaction terms in the above 
equations. 

We evaluated during the evolution of the system the mean field values, as well as the different energy 
components. Initially, the Higgs field grows towards the true vacuum and the infiaton moves towards the 
minimum of its potential and oscillates around it. We have checked that the sum of the averaged gradient, 
kinetic and potential energies (contributed by both the inflaton and the Higgs), remains constant during 
reheating, as expected, since the expansion of the universe is irrelevant in this model. We have also 
checked that the time evolution of the different energy components is the same for different lattices, 
changing the number of points N , the minimum momentum Pmin = '^i^ jT or the lattice spacing a = L/N , 
with L the lattice size. The evolution of the Higgs' v.e.v. follows three stages easily distinguished. First, 
an exponential growth of the v.e.v. towards the true vacuum. This is driven by the tachyonic instability of 
the long-wave modes of the Higgs field, that makes the spatial distribution of this field to form lumps and 
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Figure 1: The time evolution of the different types of energy (kinetic, gradient, potential, anisotropic 
components and gravitational waves for different lattices), normalized to the initial vacuum energy, after 
hybrid inflation, for a model with v — 10~^ Mp. One can clearly distinguish here three stages: tachyonic 
growth, bubble collisions and turbulence. 



bubble-like structures [231 US]. Second, the Higgs field oscillates around the true vacuum, as the Higgs' 
bubbles collide and scatter off eachother. Third, a period of turbulence is reached, during which the 
inflaton oscillates around its minimum and the Higgs sits in the true vacuum. For a detailed description 
of the dynamics of these fields see Ref. |25]. Here we will be only concerned with the details of the 
gravitational wave production. 

The initial energy density at the end of hybrid inflation is given by pq = rn^v^ /A, with iv? = Au^, so 
the fractional energy density in gravitational waves is 

'-^ = = ^l^Ah^M') > (24) 

where (^hijk^^^ , defined as a volume average like ^ J (fixhijh''^ ^ is extracted from the simulations as 

{kM')^ = y J dlogk e (Ay- i„(k)«y(t, k)uL(t, k))^^ (25) 

where Uy (t, k) is the Fourier transform of the solution of Eq. (|2ip . Then, we can compute the corre- 
sponding density parameter today (with Vl^ad — 3.5 x 10~^) 

"ow = ^ "'•f / dlogkk^ (A,,,™(k)u.,(t,k)uL(t,k)) (26) 

which has assumed that all the vacuum energy po gets converted into radiation, an approximation which 
is always valid in generic hybrid inflation models with v <C Mp, and thus iJ <C m = \f\v. 

We have shown in Fig. [T] the evolution in time of the fraction of energy density in GW. The first 
(tachyonic) stage is clearly visible, with a (logarithmic) slope twice that of the anisotropic tensor H^ . 
Then there is a small plateau corresponding to the production of GW from bubble collisions; and finally 
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Figure 2: We show here the comparison between the power spectrum of gravitational waves obtained 
with increasing lattice resolution, to prove the robustness of our method. The different realizations are 
characterized by the the minimum lattice momentum (pmin) and the lattice spacing (ma). The growth 
is shown in steps of mAt = 1 up to mt = 30, and then in and mAt = 5 steps up to mt = 60. 



there is the slow growth due to turbulence. In the next section we will describe in detail the most 
significant features appearing at each stage. 

Note that in the case that H m, the maximal production of GW occurs in less than a Hubble time, 
soon after symmetry breaking, while turbulence lasts several decades in time units of m~^. Therefore, 
we can safely ignore the dilution due to the Hubble expansion, up to times much greater than those of 
the tachyonic instability. Eventually the universe reheats and the energy in gravitational waves redshifts 
like radiation thereafter. 

To compute the power spectrum per logarithmic frequency interval in GW, fl^^{f), we just have to 
use ^TE\i. We can evaluate the power spectrum today from that obtained at reheating by converting the 
wavenumber k into frequency /. Simply using Eq. (|20p . with gr,s/go,s 100, gr,s ~ 5r and Og ~ a*, then 

/^6xlO^°Hz / ^5x 10^"HzAa1/^ (27) 

We show in Fig. [2] the power spectrum of gravitational waves as a function of (comoving) wavenumber 
k/m. We have used different lattices in order to have lattice artifacts under control, specially at late times 
and high wavenumbers. We made sure by the choice of lattice size and spacing (i.e. fc,jii,i and fc,nax) that 
all relevant scales fitted within the simulation. Note, however, that the lower bumps are lattice artifacts, 
due to the physical cutoff imposed at the initial condition, that rapidly disappear with time. We have also 
checked that the power spectrum of the scalar fields follows turbulent scaling after mt ^ 0(100), and we 
can thus estimate the subsequent evolution of the energy density distributions beyond our simulations. 



3 Lattice simulations 

The problem of determining the time evolution of a quantum field theory is an outstandingly difficult 
problem. In some cases only a few degrees of freedom are relevant or else perturbative techniques are 
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applicable. However, in our particular case, our interests are focused on processes which are necessarily 
non-linear and non-perturbative and involve many degrees of freedom. The presence of gravitational 
fields just contributes with more degrees of freedom, but does not complicate matters significantly. 

In the present paper we will use the so called classical approximation to deal with the problem. 
It consists of replacing the quantum evolution of the system by its classical evolution, for which there 
are feasible numerical methods available. The quantum nature of the problem remains in the stochastic 
character of the initial conditions. This approximation has been used with great success by several groups 
in the past [331 IMj ■ The advantage of the method is that it is fully non-linear and non-perturbative. 

Our approach is to discretize the classical equations of motion of all fields in both space and time. 
The time-like lattice spacing at must be smaller than the spatial one as for the stability of the discretized 
equations. In addition to the ultraviolet cut-off one must introduce an infrared cut-off by putting the 
system in a box with periodic boundary conditions. In this paper we have thouroughly studied a model 
with = 2A = 1/4, but we have checked that other values of the parameters do not change our results 
significantly. 

3.1 Initial conditions 

The initial conditions of the fields follow the prescription from Ref. [25] . The Higgs modes 0^ are solutions 
of the coupled evolution equations, which can be rewritten as (j>'^ + (fc^ — T)(j)k = 0, with t — M{t — tc) 
and M = {2VY^^m. The time-dependent Higgs mass follows from the initial inflaton field homogeneous 
component, Xo(^i) — Xc(l ~ Vm{ti — tc)) and Xo{ti) — —XcVm. The Higgs modes with k/M > y/n are 
set to zero, while the rest are determined by a Gaussian random field of zero mean distributed according 
to the Rayleigh distribution 

P{\<l)k\)d\cl)k\ddk = exp — — — , (28) 

V CTfc / CTfc 27r 

with a uniform random phase 9k S [0, 2tt] and dispersion given by = j/fep = P{k, Ti)/k^, where P{k, ti) 
is the power spectrum of the initial Higgs quantum fluctuations, computed in the linear approximation 
in the background of the homogeneous inflaton. In the classical limit, the conjugate momentum ^^(t) is 
uniquely determined as (/)/c(r) = F{k,T)4>k{T), with F{k,T) = lm{ifk{T)f^{T))/\fk{T)\'^, see Ref. [25] . 

The rest of the fields (the infiaton non-zero modes and the gravitational waves), are supposed to 
start from the vacuum, and therefore they are semiclassically set to zero initially in the simulations. 
Their coupling to the Higgs modes will drive their evolution, giving rise to a rapid (exponential) growth 
of the GW and infiaton modes. Their subsequent non-linear evolution will be well described by the 
lattice simulations. In the next subsections we will describe the different evolution stages found in our 
simulations. 



3.2 Tachyonic growth 

In this subsection we will compare the analytical estimates with our numerical simulations for the initial 
tachyonic growth of the Higgs modes and the subsequent growth of gravitational waves. The first check 
is that the Higgs modes grow according to Ref. [25]. There we found that 

k\Mt)\' ^v^AiT)e-''^-'^''\ (29) 

with A{t) and B{t) are given, for r > 1, as A{t) — 2r^(i/3) ' ^i''') = '^iV^ ~ 1); where 

Bi(z) is the Airy function of the second kind. We have checked that the initial growth, from mt = 6 to 
mt =10, follows precisely the analytical expression. 

The comparison between the tensor modes hij{k,t) and the numerical results is somewhat more 
complicated. We should first compute the effective anisotropic tensor Tij (k, t) from the gradients of 
the Higgs field (those of the infiaton are not relevant during the tachyonic growth), as follows, 

%(k,i) = J [V.</>'^ V,0nx,t)] , (30) 
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where Vi(/)°(x, i) = J -p^^f/s- ifZi 0°(q, e After performing the integral in x and using the delta 

function to eliminate q', we make a change of variables q ^ q + k/2, and integrate over q. Finally, with 
the use of riij(k, t), we can compute the tensor fields, 

/iy(k,i) = (167rG) [ dt' ^^"^(^~^ ) n^^.^ ^0/1^^ (k, t) = (167rG) j dt' cos k{t - t')ti,j . (31) 
Jo ^ Jo 

Using the analytic solutions one can perform the integrals and obtain expressions that agree surprisingly 
well with the numerical estimates. This allows one to compute the density in gravitational waves, Pg„, at 
least during the initial tachyonic stage in terms of analytical functions, and we reproduce the numerical 
results. We will now compare these with the analytical estimates. The tachyonic growth is dominated 
by the faster-than-exponential growth of the Higgs modes towards the true vacuum. The (traceless) 
anisotropic strees tensor 11^^ grows rapidly to a value of order ~ m?v^ , which gives a tensor 

perturbation 

\hijh'^\^'^ - 167rGw2(TOAt)2lO-3^ (32) 

and an energy density in GW, 

Pgw/Po M-kGv'^ {mMflQ-^ - Gw^ , (33) 

for mAi ~ 16. In the case at hand, with v — lO^'^Afp, we find p^^j/Po ^ 10^^ at symmetry breaking, 
which coincides with the numerical simulations at that time, see Fig. [T] 

As shown in Ref. [25], the spinodal instabilities grow following the statistics of a Gaussian random 
field, and therefore one can use the formalism of [41j to estimate the number of peaks or lumps in the 
Higgs spatial distribution just before symmetry breaking. As we will discuss in the next section, these 
lumps will give rise via non-linear growth to lump invagination and the formation of bubble- like structures 
with large density gradients, expanding at the relativistic speeds and colliding among themselves giving 
rise to a large GWB. The size of the bubbles upon collision is essentially determined by the distance 
between peaks at the time of symmetry breaking, but this can be computed directly from the analysis of 
Gaussian random fields, as performed in Ref. 25J. This analysis works only for the initial (linear) stage 
before symmetry breaking. Nevertheless, we expect the results to extrapolate to later times since once a 
bubble is formed around a peak, it remains there at a fixed distance from other bubbles. This will give 
us an idea of the size of the bubbles at the time of collision. 



3.3 Bubble collisions 

The production of gravitational waves in the next stage proceeds through 'bubble' collisions. In Ref. [24] 
we showed that during the symetry breaking, the Higgs field develops lumps whose peaks grow up to a 
maximum value |0|max/i' = 4/3, and then decrease creating approximately spherically symmetric bubbles, 
with ridges that remain above |0| = v. Finally, neighboring bubbles collide and high momentum modes 
are induced via field inhomogeneities. Since initially only the Higgs field sources the anisotropic stress- 
tensor Hij, then we expect the formation of structures in the spatial distribution of the GW energy 
density correlated with the Higgs lumps. In this sub-section we will give an estimate of the burst in GW 
produced by the first collisions of the Higgs bubble-like structures. 

As for the collision of vaccum bubbles in first order phase transitions [S] , we can give a simple estimate 
of the order of magnitude of the energy fraction radiated in the form of gravitational waves when two 
Higgs bubble- like structures collide. A similar stimation is indeed presented in [151 [S^. In general, the 
problem of two colliding bubbles has several time and length scales: the duration of the collision. At; the 
bubbles' radius R at the moment of the collision; and the relative speed of the bubble walls. The typical 
size of bubbles upon collisions, is of the order of i? « 10m~^, while the growth of the bubble's wall is 
relativistic, see Ref. [25]. Then we can assume than the time scale associated with bubble collisions is 
also At ~ R. Assuming the bubble walls contain most of the energy density, it is expected that the 
asymmetric collisions will copiously produce GW. 

Far from a source that produces gravitational radiation, the dominat contribution to the amplitude 
of GW is given by the acceleration of the quadrupole moment of the Higgs field distribution. Given the 



10 



energy density of the Higgs field, pn, we can compute the (reduced) quadrupole moment of the Higgs 
field spatial distribution, Qij = J dPx (xiXj — x^Sij/i) pnix), such that the amplitude of the gravitational 
radiation, in the TT gauge, is given by hij ~ {2G/r)Qij . A significant amount of energy can be emitted in 
the form of gravitational radiation whenever the quadrupole moment changes significantly fast: through 
the bubble collisions in this case. The power carried by these waves can be obtained via ([T7]) as 

Po^ = ^ j d^l {q,,W) . (34) 

Omitting indices for simplicity, as the power emitted in gravitational waves in the quadrupole approxi- 
mation is of order ^ G{Q)^, while the quadrupole moment is of order Q ^ R^pn, we can estimate 
the power emitted in GW upon the collision of two Higgs bubbles as 

Pow -gI^^^' ^ GplR' (35) 

The fraction of energy density carried by these waves, p^^ ^ P^^^At/R^ ~ P^^/R'^ ^ Gp^R^, 
compared to that of the initial energy stored in the two bubble-like structures of the Higgs field, will be 
Pgw/ph = GpuR^. Since the expansion of the universe is negligible during the bubble collision stage, the 
energy that drives inflation, po ~ m'^v'^, is transferred essentially to the Higgs modes during preheating, 
within an order of magnitude, see Fig. [T] Thus, recalling that R ~ 10m~^, the total fraction of energy 
in GW produced during the bubble collisions to that stored in the Higgs lumps formed at symmetry 
breaking, is given by 

O.lGpoR^ - (v/Mpf, (36) 



Po 



giving an amplitude which is of the same order as is observed in the numerical simulations, see Fig.[T] Of 
course, an exhaustive analytical treatment of the production of GW during this stage of bubble collisions 
remains to be done, but we leave it for a future publication. 



3.4 Turbulence 

The development of a turbulent stage is expected from the point of view of classical fields, as turbulence 
usually appears whenever there exists an active (stationary) source of energy localized at some scale fcjn 
in Fourier space. The oscillating infiaton zero-mode plays the role of the pumping-energy source, acting 
at a well defined scale fcin in Fourier space, given by the frequency of the infiaton oscillations. Apart from 
k\m there is no other scale in Fourier space where energy is accummulated, dissipated and/or infused. So, 
as turbulence is characterized by the transport of some conserved quantity, energy in our case, we should 
expect a fiow of energy from km towards higher (direct cascade) or smaller (inverse cascade) momentum 
modes. In typical turbulent regimes of classical fluids, there exits a sink in Fourier space, corresponding to 
that scale at which the (direct) cascade stops and energy gets dissipated. However, in our problem there 
is no such sink so that the transported energy cannot be dissipated, but instead it is used to populate 
high-momentum modes. For the problem at hand, there exists a natural initial cut-off /cout ^ A^/^w, 
such that only long wave modes within k < A:out7 develop the spinodal instability. Eventually, after the 
tachyonic growth has ended and the flrst Higgs' bubble- like structures have collided, the turbulent regime 
is established. Then the energy flows from small to greater scales in Fourier space, which translates into 
the increase of /cout in time. 

When the turbulence has been fully established, if the wave (kinetic) turbulence regime of the fields' 
dynamics is valid, the time evolution of the variance of a turbulent field /(x, i), should follow a power- 
law-like scaling [13] 

Var(/(i)) = (/(t)2) - (/(t)) V t-'P, (37) 

with p — 1/{2N — 1) and the number of scattering fields in a 'point-like collision'. In Fig. [3] we have 
plotted the time evolution of the variances of the Infiaton x and of the Higgs modulus cj) — \/J2a '^a' ''^^'^ 
fitted the data with a power-law like p7p , obtaining 
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Figure 3: Variance of the Inflaton and the Higgs field as a function of time, the former normalized to 
its critical value, the latter normalized to its v.e.v.. As expected in a turbulent regime, these variances 
follow a power law ^ t~'^P withp a certain critical exponent, although the slope of the Inflaton's variances 
evolves in time. The curves are produced from an average over 10 different statistical realizations. 
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where the last brackets on the right correspond to the range in time (in units of m^^) for which we 
fitted the data. As can be seen in Fig.[3l the slope of the Higgs field (in logarithmic scale), 2p^ ^ 2/7, 
remains approximately constant in time, corresponding to a 4-ficld dominant interaction. However, the 
slope of the Inflaton's variance increases in time, i.e. the critical exponent Pj of the Inflaton decreases, 
until it reaches a stationary stage at mt 100. We will not try to explain here the origin of such an 
effective critical exponents as extracted from the simulations. We will just stress that we have checked 
the robustness of those values under different lattice configurations (A^, Pmin) and different statistical 
realizations. Actually, when turbulence has fully developed, it is expected that the distribution function 
of the classical turbulent fields, the inflaton and the Higgs here, follow a self-similar evolution [43] 



(38) 



with p the critical exponent of the fields' variances and 7 a certain factor ~ 0(1), which depends on 
the type of turbulence developed. Looking at ([551) . we see that the exponent p determines the speed 
of the particles' distribution in momentum space: given a specific scale kc that scale evolves in time 
as fcc(t) = kc{to){t/to)P . In the simulations, we have seen that the evolution of the Higgs occupation 
number follows Eq. ((38|) with p « 1/7, as expected from the Higgs variance, and 7 « 2.7. Whereas the 
evolution of the Inflaton occupation number follows (I38p even more accurately than the Higgs, with an 
"effective" exponent (once the asymptotic regime is achieved) p « 1/5, and 7 ~ 3.9. In Figs. [4] we have 
plotted the occupation numbers of the Higgs and the Inflaton, also inverting the relation of Eq. ([55)1 in 
order to extract the universal time- independent no{k) functions of each field. As shown in those flgures, 
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Figure 4: Some snapshots of the evolution of the spectral particle occupation numbers of the Higgs and 
the Inflaton fields at different times, each averaged over 10 statistical realizations. We multiply them by 
k'^ so we can see better the scaling behaviour. In the upper right corner, we plot the inverse relation 
of no{kt~P) = t'^Pn{k,t), also averaged over 10 realizations for each time. The scaling behaviour 
predicted by wave kinetic turbulent theory [43] , is clearly verified. 



the distributions follow nicely the expected scaling behaviour. The universal functions nQ^k) plotted in 
Figs, m have been obtained from averaging over ten statistical realizations for each time. 

The advantage of the development of a turbulence behaviour is obvious: it allows us to extrapolate 
the time evolution of the fields' distributions till later times beyond the one we can reach with the 
simulations. Moreover, the fact that the turbulence develops so early after the tachyonic instability, also 
allow us to check for a long time of the simulation, the goodness of the description of the dynamics of the 
fields, given by the turbulent kinetic theory developed in Ref. [43]. We have fitted the averaged universal 
functions no{k) with expressions of the form fc*no(fc) = P{k)e~'^^'^\ with P{k) and Q{k) polynomials in 
k. There is no fundamental meaning associated with such a fit, but it is very useful to have an analytical 
control over no(fc), since this allows us to track the time-evolution of n{k, t) through Eq. ([38]) . Actually, 
the classical regime of the evolution of some bosonic fields ends when the system can be relaxed to the 
Bose-Einstein distribution. Since we cannot reach that moment, we can at least estimate the moment in 
which the initial energy density gets fully transferred to the Higgs classical modes. Using Ea. ([551) and 
the fit to the universal rio(fc) of the Higgs, we find that the initial energy density is totally transfered to 
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Figure 5: Time evolution of the GW spectra from mt = 6 to mt = 2000. The amphtude of the spectra 
seems to saturate after mt ~ 100, although the high momentum tail still moves slowly to higher values 
of k during the turbulent stage. 

the Higgs when (in units rn = 1) 



(39) 



where we have assumed that the Higgs' modes have energy Ek{k, t) — k n{k, t). In our case, with A = 1/8, 
the conversion of the initial energy density into Higgs particles and therefore into radiation is complete 
by t ~ 6 X 10'*m~^. Therefore, if we consider this value as a lower bound for the time that classical 
turbulence requires to end, we see that turbulence last for a very long time compared to the time-scale 
of the initial tachyonic and bubbly stages. Thus, if GW were significatively sourced during turbulence, 
one should take into account corrections from the expansion of the universe. 

In Fig. [HI we show the evolution of the GW spectra up to times mt = 2000, for a lattice of (N,pmin) 
— (128,0.15). It is clear from that figure that the amplitude of the GW saturates to a value of order 
Pgw/Po ~ 2 ■ 10~^. At mt w 50, the maximum amplitude of the spectra has already reached Pqw/Po ~ 
10~^, while at time mt « 100, the maximum has only grown a factor of 2 with respect to mt « 50. 
From times mt « 150 till the maximum time we reached in the simulations, mt = 2000, the maximum of 
the amplitude of the spectrum does not seem to change significantly, slowly increasing from « 2 • 10~^ 
to « 2.5 • 10"^. Despite this saturation, we see in the simulations that the long momentum tail of the 
spectrum keeps moving towards greater values. This displacement is precisely what one would expect from 
turbulence, although it is clear that the amplitude of the new high momentum modes never exceed that of 
lower momentum. In order to disscard that this displacement towards the UV is not a numerical artefact, 
one should further investigate the role played by the turbulent scalar fields as a source of GW. Here, we 
just want to remark that the turbulent motions of the scalar fields, seem not to increase significatively 
anymore the total amplitude of the GW spectrum. Indeed, in a recent paper [22] where GW production 
at reheating is also considered, it is stated that GW production from turbulent motion of classical scalar 
fields, should be very supressed. That is apparently what we observe in our simulations although, as 
pointed above, this issue should be investigated in a more detailed way. Anyway, here we can conclude 
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that the expansion of the Universe during reheating in these hybrid models, does not play an important 
role during the time of GW production, and therefore we can be safely ignore it. 



4 Gravitational Waves from Chaotic Inflation 

The production of a relic GWB at reheating was first addressed by Khlebnikov and Tkachev (KT) in 
Ref. [18j . both for the quadratic and quartic chaotic inflation scenarios. Recently, chaotic scenarios were 
revisited in Ref. [20l [21] . Also very recently, Ref. [22] studied in a very detail way, the evolution of GW 
produced at preheating in the case of a massless inflaton with an extra scalar field. 

In Refs. [T^ and [501, the procedure to compute the GW from reheating relied on Weinberg's formula 
for flat space-time [45]. However, in chaotic models, the expansion of the universe cannot be neglected 
during reheating, so Weinberg's formula can only be used in an approximated way, if the evolution of 
the universe is considered as an adiabatic sequence of stationary universes. In Ref. [17] . however, we 
adopted a different approach that takes into account the expansion of the universe in a self-consistent 
manner, and allows us to calculate at any time the energy density and power spectra of the GW produced 
at reheating (see section 2). Using our technique, we will show in this section that we reproduce, for 
specific chaotic models, similar results to those of other authors. In particular, we adapted the publicly 
available LATTICEEASY code [31], taking advantage of the structure of the code itself, incorparating 
the evolution of Eq. ([7|), together with the equations of the scalar fields, Eqs. ((4|), into the staggered 
leapfrog integrator routine. 

Here we will concentrate only in an scenario with a massless inflaton x, either accompanied or not by 
an extra scalar field 4>. Such scenarios are described by the potential 

- ^X' + (40) 
Rescaling the time by and the physical fields by a conformal transformation as 

dr^^-^.{Q)V\dt, (41) 
a(0)x(0) a(0)x(0) a{Q)^^ ' ' ^ ' 

then the equations of motion of the inflaton and of the extra scalar field, Eq. Q, can be rewritten in 
terms of the conformal variables as 



a" 



Xc - ^ xc - —xc + ixi + qrc)xc = (42) 

a 

c^'^-V^,~—Xc + qxl<Pc = 0, (43) 
a 

where the prime denotes derivative with respect to conformal time. Since the universe expands as 
radiation-like in these scenarios, a{T) ^ r, so the terms proportional to a" /a in Eqs. (I42p and (|43|) 
are soon negligible, as explicitly checked in the simulations. Thanks to this, the model is conformal to 
Minkowski. 

The parameter q = g^/X controls the strength and width of the resonance. For the case of a massless 
inflaton without an extra scalar field, we just set g = in Eq. ([42]) and ignore Eq. (|43)) . However, in that 
case, fiuctuations of the infiaton also grow via parametric resonance. Actually, they grow as if they were 
fiuctuations of a scalar field coupled to the zero- mode of the inflaton with effective couplig q = /X = 3, 
see Ref. [35]. Following Refs. [H] and |2U], we set A = 10~^^ and q = 120. Since this case is also computed 
in [22], we can also compare our results with theirs. Moreover, we also present results for the pure Ax^ 
model with no extra scalar field, a case only shown in Ref. [TF. 

We begin our simulations at the end of infiation, when the homogeneous infiaton verifies xo ~ 0.342Mp 
and Xo ~ 0. We took initial quantum (conformal) fiuctuations 1/V2k for all the modes up to a certain 
cut-off, and only added an initial zero-mode for the inflaton, Xc(0) = 1, Xc(0)' = 0. In Figs. [6] and [T] we 
show the evolution of fi^^ during reheating, normalized to the instant density at each time step, for the 
coupled and the pure case, respectively. In the case with an extra scalar field, the amplitude of the GWB 
saturates at the end of parametric resonance, when the fields variances have been stabilized. This is the 
beginnig of the turbulent stage in the scalar fields, which seems not to source anymore the production of 
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Figure 6: The spectrum of the gravitational waves' energy density, for coupled case with A = 10^^^ and 
g^/A = 120. The spectrum is shown accumulated up to different times during GW production, so one 
can see its evolution. At each time, it is normalized to the total instant density. This plot corresponds 
to a N = 128 lattice simulation, from t = to r = 240. 



GWs, as already stated in Refs.|20'i'22]. For the pure case, we also see the saturation of the amplitude of 
the spectra, see Fig. [71 although the high momentum tail seems to slightly move toward higher values. 

Of course, in either case, with and without an extra field (fi, in order to predict today's spectral window 
of the GW spectrum, we have first to normalize their energy density at the end of GW production to the 
total energy density at that moment; then to redshift the GW spectra from that moment of reheating, 
taking into account that the rate of expansion have changed significantly since the end of inflation, see 
Eg. ([20]) . In particular, the shape and amplitude of GW spectra for the case with the extra scalar field 
coupled to the inflaton with q — 120, seems to coincide with the espectra shown in Ref. [22]. On the 
other hand, we also reproduce a similar spectra to the one shown in TSj, for the case of the pure quartic 
model. Thanks to the tremendous gain in computer power, we were able to resolve the 'spiky' pattern 
of that spectrum with great resolution. For the first time, it is clearly observed the exponential tail for 
large frequencies, not shown in Ref. [18]. The most remarkable fact, is that we also confirm that the 
peak structure in the GW power spectrum, see Fig. [71 remains clearly visible at times much later than 
the one at which those peaks have dissapeared in the scalar fields' power spectrum. So, as pointed out 
in Ref. [18], this characteristic feature will allow us to distinguish this particular model from any other. 

Let us emphasize that we have run the simulations till times much greater than that of the end of 
the resonance stage, both for the pure and the coupled case. The role of the turbulence period after 
preheating seems, therefore, not to be very important, despite its long duration. Apparently, the no-go 
theorem about the suppresion of GW at turbulence, discussed in [55], is fulfilled. In Refs. [571 113 it was 
pointed out that gauge couplings or trilinear interactions could be responsible for a fast thermalization 
of the universe after infiation (see also Ref. [49]), but as long as this takes place after the end of the 
resonace stage, in principle this should not affect the results shown above. 



5 Conclusions 



To summarize, we have shown that hybrid models are very efficient generators of gravitational waves at 
preheating, in three well defined stages, first via the tachyonic growth of Higgs modes, whose gradients 
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Figure 7: The spectrum of the gravitational waves' energy density, for the pure case, with A = 10^^"'. 
Again, we show the spectrum accumulated up to different times during GW production, normalized to 
the total instant density at each time. The plot corresponds to a N = 128 lattice simulation, from t = 



act as sources of gravity waves; then via the collisions of highly relativistic bubble-like structures with 
large amounts of energy density, and finally via the turbulent regime (although this effect does not seem 
to be very significant in the presence of scalar sources) , which drives the system towards thermalization. 
These waves remain decoupled since the moment of their production, and thus the predicted amplitude 
and shape of the gravitational wave spectrum today can be used as a probe of the reheating period in 
the very early universe. The characteristic spectrum can be used to distinguish between this stochastic 
background and others, like those arising from NS-NS and BH-BH coalescence, which are decreasing with 
frequency, or those arising from inflation, that are flat [50] . 

We have plotted in Fig. [8] the sensitivity of planned GW interferometers like LIGO, LISA and BBO, 
together with the present bounds from CMB anisotropies (GUT inflation) , from Big Bang Nucleosynthesis 
(BBN) and from milisecond pulsars (ms pulsar). Also shown are the expected stochastic backgrounds of 
chaotic inflation models like At/i*, both coupled and pure, as well as the predicted background from two 
different hybrid inflation models, a high-scale model, with v — lQ~^Mp and X ^ ^ 0.1, and a low-scale 
model, with v = 10~^Mp and A ~ ~ 10~^^, corresponding to a rate of expansion H ~ 100 GeV. The 
high-scale hybrid model produces typically as much gravitational waves from preheating as the chaotic 
inflation models. The advantage of low-scale hybrid models of inflation is that the background produced 
is within reach of future GW detectors like BBO [4]. It is speculated that future high frequency laser 
interferometers could be sensitive to a GWB in the MHz region [T^, although they are still far from the 
bound marked with an interrogation sign. 

For a high-scale model of inflation, we may never see the predicted GW background coming from pre- 
heating, in spite of its large amplitude, because it appears at very high frequencies, where no detector has 
yet shown to be sufficiently sensitive, unless the spectrum can be extrapolated to lower frequencies, where 
there are interferometric detectors like BBO which could see a signal. On the other hand, if inflation oc- 
cured at low scales, even though we will never have a chance to detect the GW produced during inflation in 
the polarization anisotropies of the CMB, we do expect gravitational waves from preheating to contribute 
with an important background in sensitive detectors like BBO. The detection and characterization of such 
a GW background, coming from the complicated and mostly unknown epoch of rehating of the universe, 



to r = 2000. 
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Figure 8: The sensitivity of the different gravitational wave experiments, present and future, compared 
with the possible stochastic backgrounds; we include the White Dwarf Binaries (WDB) (4^ and chaotic 
preheating (A0^, coupled and pure) for comparison. Note the two well differentiated backgrounds from 
high-scale and low-scale hybrid inflation. The bound marked (?) is estimated from ultra high frequency 
laser interferometers' expectations [T2] . 

may open a new window into the very early universe, while providing a new test on inflationary cosmology. 



Acknowledgments 

We wish to thank Alfonso Sastre for a very fruitful collaboration. JGB thanks the organisers of the 
Japan GRG 17th Meeting in Nagoya for a very enjoyable conference and legendary hospitality. This 
work is supported in part by CICYT projects FPA2003-03801 and FPA2006-05807, by EU network 
"UniverseNet" MRTN-CT-2006-035863 and by CAM project HEPHACOS S-0505/ESP-0346. D.G.F. 
and A.S. acknowledges support from a FPU-Fellowship from the Spanish M.E.C. We also acknowledge 
use of the MareNostrum Supercomputer under project AECT-2007-1-0005. 



References 

[1] M. Maggiore, Phys. Rept. 331, 283 (2000); C. J. Hogan, |arXiv:astro-ph/0608567| A. Buonanno, 
"Gravitational waves," larXiv:0709."4682l [gr-qc]. 



[2] B. Abbott et al. [LIGO Scientific Collaboration], Astrophys. J. 659, 
[arXiv:astro-ph/0608606 . LIGO Home Page: |http : //www . ligo ■ caltech . edu/| 



918 (2007) 



[3] S. A. Hughes. !arXiv:0711.0188l [gr-qc] . LISA Home Page: |http://lisa.esa.int 



[4] V. Corbin and N. J. Cornish, Class. Quant. Grav. 23, 2435 (2006); G. M. Harry, P. Fritschel, 
D. A. Shaddock, W. Folkncr and E. S. Phinney, Class. Quant. Grav. 23, 4887 (2006). BBO Home 
Page: ,http : / / universe . nasa . gov/ new/ progr am/bbo . html, 



18 



S. Kawamura et ai, Class. Quant. Grav. 23, S125 (2006). 

T. L. Smith, E. Pierpaoli and M. Kamionkowski, Phys. Rev. Lett. 97, 021301 (2006). 



B. Allen, gr-qc/9604033 



I. H. Stairs, Living Rev. Rel. 6, 5 (2003). 

A. Kosowsky, M. S. Turner and R. Watkins, Phys. Rev. Lett. 69, 2026 (1992), Phys. Rev. D 45, 
4514 (1992); A. Kosowsky and M. S. Turner, Phys. Rev. D 47, 4372 (1993); M. Kamionkowski, 
A. Kosowsky and M. S. Turner, Phys. Rev. D 49, 2837 (1994). 

A. Nicolis, Class. Quant. Grav. 21, L27 (2004); C. Grojean and G. Servant, "Gravitational waves 
from phase transitions at the electroweak scale and beyond," arXiv:hep-ph/0607107, 

A. Kosowsky, A. Mack and T. Kahniashvih, Phys. Rev. D 66, 024030 (2002); A. D. Dolgov, D. Grasso 
and A. Nicohs, Phys. Rev. D 66, 103505 (2002). G. Gogoberidze, T. Kahniashvih and A. Kosowsky, 
"The spectrum of gravitational radiation from primordial turbulence," arXiv:0705.1733l [astro-ph]. 

A. Nishizawa et ai, "Laser-interferometric Detectors for Gravitational Wave Background at 100 
MHz : Detector Design and Sensitivity," arXiv:0710.1944 [gr-qc]; A. M. Cruise and R. M. J. Ingley, 
Class. Quant. Grav. 23, 6185 (2006). 

F. Pegoraro, L. A. Radicati, P. Bernard and E. Picasso, Phys. Lett. A 68, 165 (1978); G. Gcmme, 



A. Chincarini, R. Parodi, P. Bernard and E. Picasso, gr-qc/0112021; R. Ballantini et ai, 
g r-qc/0502054 , 

A. A. Starobinsky, JETP Lett. 30, 682 (1979) [Pisma Zh. Eksp. Tcor. Fiz. 30, 719 (1979)]. 

L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. Lett. 73, 3195 (1994); L. Kofman, 
A. D. Linde and A. A. Starobinsky, Phys. Rev. D 56, 3258 (1997). 

J. Garcia-Belhdo and D. G. Figueroa, Phys. Rev. Lett. 98, 061302 (2007) [arXiv:astro-ph/070101^ . 

J. Garcia-Bellido, D. G. Figueroa and A. Sastre, "A Gravitational Wave Background from Reheating 
after Hybrid Inflation," larXiv:0707."0 839 [hep-ph], to appear in PRD. 

S. Y. Khlebnikov and I. I. Tkachev, Phys. Rev. D 56, 653 (1997). 

J. Garcia-Bellido, "Preheating the universe in hybrid inflation," Proceedings of the 33rd Ren- 
contres de Moriond: Fundamental Parameters in Cosmology, Les Arcs, France, 17-24 Jan 1998, 
|arXiv:hep-ph/ 9804205, 

R. Easther and E. A. Lim, JCAP 0604, 010 (2006). 

R. Easther, J. T. . Giblin and E. A. Lim, "Gravitational Wave Production At The End Of Inflation," 
arXiv :astro-ph/061229 4, 

J. F. Dufaux, A. Bergman, G. N. Felder, L. Kofman and J. P. Uzan. 'a TXiv:0707.0875 ' [astro-ph]. 

A. D. Linde, Phys. Rev. D 49, 748 (1994); J. Garcia-Bellido and A. D. Linde, Phys. Rev. D 57, 6075 

(1998) . 

G. N. Felder, J. Garcia-Bellido, P. B. Greene, L. Kofman, A. D. Linde and I. Tkachev, Phys. Rev. 
Lett. 87, 011601 (2001); G. N. Felder, L. Kofman and A. D. Linde, Phys. Rev. D 64, 123517 (2001). 

J. Garcia-Belhdo, M. Garcia Perez and A. Gonzalez- Arroyo, Phys. Rev. D 67, 103501 (2003); 

J. Garcia-Belhdo and E. Ruiz Morales, Phys. Lett. B 536, 193 (2002). 

J. Garcia-Belhdo, D. Y. Grigoriev, A. Kusenko and M. E. Shaposhnikov, Phys. Rev. D 60, 123504 

(1999) . J. Garcia-Belhdo, M. Garcia Perez and A. Gonzalez- Arroyo, Phys. Rev. D 69, 023504 (2004); 
A. Tranberg and J. Smit, JHEP 0311, 016 (2003). 



19 



[28] A. Diaz-Gil, J. Garcia-Bellido, M. Garcia Perez and A. Gonzalez- Arroyo, PoS LAT2005, 242 (2006); 
PoS LAT2007, 052 (2007), arXiv:0710.0580 [hep-lat]; "Magnetic field production during preheating 
at the electroweak scale," arXiv:0712.4263 [hep-ph]. 

[29] MareNostrum Supercomputer Home Page: 
|http : //www ■ bsc ■ es/ 

[30] IFT cluster Home Page: 

|http : //lattice . f t . uam . es/if tcluster/| 

[31] G. N. Felder and I. Tkachev, arXiv: hep-ph/0011159! 
http : //www . science . smith . edu/ depar tments/, 
Physics/f staff /gfelder/latticeeasy/ 

[32] J. Garcia-BeUido, A. D. Linde and D. Wands, Phys. Rev. D 54, 6040 (1996). 

[33] R. AUahverdi, K. Enqvist, J. Garcia-BeUido and A. Mazumdar, Phys. Rev. Lett. 97, 191304 (2006); 
R. Allahverdi, K. Enqvist, J. Garcia-Bellido, A. Jokinen and A. Mazumdar, JCAP 0706 (2007) 019 
| arXiv:heFih70610134 . 

[34] S. Y. Khlebnikov and I. I. Tkachev, Phys. Rev. Lett. 77, 219 (1996); 79, 1607 (1997); T. Prokopec 
and T. G. Roos, Phys. Rev. D 55, 3768 (1997). 

[35] V. F. Mukhanov, H. A. Feldman and R. H. Brandenberger, Phys. Rept. 215, 203 (1992). 

[36] S. Carroll, ^^Spacetime and Geometry: An introduction to General Relativity" Addison Wesley 
(2003). 

[37] K. G. Wilson, Phys. Rev. D 10, 2445 (1974). 

[38] F. Coo per, S. Habi b, Y. Kluger, E. Mottola, J. P. Paz and P. R. Anderson, Phys. Rev. D 50, 2848 
(1994) |arXiv:hep-ph/9405352|; F. C ooper, S. Habib, Y. Kluger and E. Mottola, Phys. Rev. D 55, 
6471 (1997) [a rXiv:hep-ph/961034^ . 

[39] D. Boyanovsky, D. Cormier, H. J. de Vega, R. Holman, A. Singh and M. Srednicki, Phys. Rev. D 
56, 1939 (1997) [arXiv:hep-ph/9 703327]; D. Boyanovsky , H. J. de Vega, R. Holman and J. Salgado, 
Phys. Rev. D 59, 125009 (1999)^ |arXiv:hep-ph/9811273| . 

[40] J. Baacke, K. Heitmann and C. Pat zold, Phys. Rev. D 55 , 2320 (1997) |larXiv:hep-th/9608006] ; 
Phys. Rev. D 55, 7815 (1997) |arXiv:hep-ph/9612264| ; Phys. Rev. D 56, 6556 (1997) 
[arXiv:hep-ph/970 6274 . 

[41] J. M. Bardeen, J. R. Bond, N. Kaiser and A. S. Szalay, Astrophys. J. 304, 15 (1986). 

[42] G. N. Felder and L. Kofman, Phys. Rev. D 75 (2007) 043518 [ arXiv:hep-ph/0606256 '. 

[43] R. Micha and L L Tkachev, Phys. Rev. Lett. 90, 121301 (2003); Phys. Rev. D 70, 043538 (2004). 

[44] |http://lattice.ft.uam.es/gw 

[45] S. Weinberg, "Gravitation and Cosmology" , John Wiley & Sons, 1972. 

[46] P. B. Greene, L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. D 56, 6175 (1997) 
|arXiv:hep-ph/9705347, . 

[47] A. J. Farmer and E. S. Phinney, Mon. Not. Roy. Astron. Soc. 346, 1197 (2003). 

[48] J. F. Dufaux, G. N. Felder, L. Kofman, M. Peloso and D. Podolsky, JCAP 0607, 006 (2006). 

[49] G. N. Felder and L. Kofman, Phys. Rev. D 63, 103503 (2001). 

[50] T. L. Smith, M. Kamionkowski and A. Cooray, Phys. Rev. D 73, 023504 (2006). 



20 



